BARRETO et al. 2022: Habitat loss and biodiversity gain (BIOCON-S-22-00746)
Model formulas
For each of our response variables we run three candidate models for each scale:
Model 0 <- ~1, absence of effect
Model 1 <- percentage of forest cover at 3 and 5km radius
Model1.2 <- non-linear forest cover at 3 and 5km radius (quadratic model)
Loading data and gathering into one datasets calculated in script ‘data_wrang_GIT.Rmd’, content: 1) “beetles”: raw data for species occurrence per landscape/site; 2) “hab.data”: data containing classes of habitat association; 3) “env.data”: environmental variables, forest cover percent at 3 and 5km radius buffer 4) “div_hab”: diversity (alpha, beta and gamma) between classification habitat association 5) “rc.df”: Beta RC calculations 6) “ab_hab”: mean and median abundance between groups of habitat association 7) “st.hab_ab”: abundance data between classification habitat association per site
| Landscape | gamma_FS | gamma_NFS | alpha_FS | alpha_NFS | betad_FS | betad_NFS | ab.land_FS | ab.land_NFS | brc.fsNA | brc.nfsNA | brc.abund.fsNA | brc.abund.nfsNA | fc_3km | fc_5km |
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| 148 | 11 | 9 | 3.75 | 2.25 | 7.25 | 6.75 | 66 | 40 | -0.2917857 | 0.1800000 | 0.0007150 | 0.4310978 | 38.52523 | 43.86337 |
| 215 | 19 | 7 | 5.62 | 2.00 | 13.38 | 5.00 | 174 | 51 | 0.1050000 | -0.3104762 | 0.2842485 | 0.3969207 | 30.15485 | 25.23160 |
| 263 | 16 | 13 | 6.00 | 4.50 | 10.00 | 8.50 | 393 | 214 | -0.4914286 | -0.0328571 | 0.6290100 | 0.3908194 | 10.10641 | 13.41180 |
| 266 | 11 | 8 | 3.25 | 1.38 | 7.75 | 6.62 | 70 | 12 | -0.2290476 | 0.1309524 | 0.1069164 | 0.1522952 | 30.87525 | 31.16090 |
| 282 | 17 | 11 | 6.50 | 4.12 | 10.50 | 6.88 | 219 | 212 | -0.3271429 | -0.1109524 | 0.2140712 | 0.7270127 | 15.26270 | 16.35308 |
| 291 | 9 | 3 | 3.71 | 1.71 | 5.29 | 1.29 | 150 | 49 | -0.3876190 | -0.5861905 | 0.0977168 | -0.4293818 | 45.74294 | 45.13213 |
| 317 | 12 | 4 | 5.00 | 1.38 | 7.00 | 2.62 | 152 | 56 | -0.3832143 | -0.5346429 | 0.1839339 | -0.4355784 | 48.75805 | 47.85138 |
| 329 | 21 | 14 | 8.62 | 5.75 | 12.38 | 8.25 | 369 | 262 | -0.0821429 | 0.2232143 | 0.9174174 | 0.4880952 | 14.73791 | 13.62588 |
| 333 | 14 | 15 | 4.50 | 4.62 | 9.50 | 10.38 | 158 | 219 | -0.1314286 | 0.1921429 | 0.6792864 | 0.6124339 | 13.71403 | 14.50960 |
| 335 | 15 | 8 | 5.75 | 2.38 | 9.25 | 5.62 | 140 | 113 | -0.2553571 | -0.3357143 | 0.3477763 | -0.9467563 | 23.87913 | 19.21670 |
| 359 | 16 | 17 | 6.00 | 5.50 | 10.00 | 11.50 | 149 | 133 | -0.2253571 | 0.3067857 | 0.7706278 | 0.8152081 | 18.10208 | 17.62431 |
| 399 | 20 | 14 | 7.88 | 5.12 | 12.12 | 8.88 | 432 | 169 | -0.1242857 | 0.2692857 | 0.9614376 | 0.7002717 | 29.73616 | 24.97851 |
RESULTS
We found 28 species to be forest specialists, those that come from the Atlantic Forest domain and exclusively collected in forested habitats; while 23 were non-forest specialists, species from a broader domain and/or less restrictive to forest habitats (hereafter FS and NFS, respectively).
FOREST SPECIALISTS
Gamma diversity
Modelling
g0 <- glm(gamma_FS ~ 1, data=data_div, family= Gamma (link = "identity"))
g5k <- glm(gamma_FS ~ fc_5km, data=data_div, family= Gamma (link = "identity"))
g5kq <- glm(gamma_FS ~ fc_5km + I(fc_5km^2), data=data_div, family= Gamma (link = "identity"))
g3k <- glm(gamma_FS ~ fc_3km, data=data_div, family= Gamma (link = "identity"))
g3kq <- glm(gamma_FS ~ fc_3km + I(fc_3km^2), data=data_div, family= Gamma (link = "identity"))
## AIC-based model selection
(g.fs.aic <- AICctab(g0, g5k, g5kq, g3k, g3kq,
base=T,weights=T, logLik=T))## logLik AICc dLogLik dAICc df weight
## g5k -27.9 64.8 4.7 0.0 3 0.731
## g3k -29.6 68.2 3.0 3.4 3 0.136
## g5kq -27.8 69.4 4.8 4.6 4 0.074
## g0 -32.6 70.6 0.0 5.8 2 0.041
## g3kq -29.2 72.2 3.4 7.3 4 0.019
Summary top model
summary(g5k)##
## Call:
## glm(formula = gamma_FS ~ fc_5km, family = Gamma(link = "identity"),
## data = data_div)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -0.23181 -0.12743 -0.03948 0.15159 0.27856
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 20.68384 1.96964 10.50 1.01e-06 ***
## fc_5km -0.21451 0.05782 -3.71 0.00404 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for Gamma family taken to be 0.0352002)
##
## Null deviance: 0.74994 on 11 degrees of freedom
## Residual deviance: 0.34328 on 10 degrees of freedom
## AIC: 61.839
##
## Number of Fisher Scoring iterations: 3
g.fs <- g5k# rename top model to save and plot belowInspect residuals:
sim <- simulateResiduals(g5k)
plot(sim)testDispersion(sim)##
## DHARMa nonparametric dispersion test via sd of residuals fitted vs.
## simulated
##
## data: simulationOutput
## dispersion = 1.1487, p-value = 0.632
## alternative hypothesis: two.sided
hnp::hnp(g5k)## Gamma model
boot::glm.diag.plots(g5k) ## for glm poisson or neg binomPrediction:
Alpha diversity
Modelling
a0 <- glm(alpha_FS ~ 1, data=data_div, family= Gamma (link = "identity"))
a5k <- glm(alpha_FS ~ fc_5km, data=data_div, family= Gamma (link = "identity"))
a5kq <- glm(alpha_FS ~ fc_5km + I(fc_5km^2), data=data_div, family= Gamma (link = "identity"))
a3k <- glm(alpha_FS ~ fc_3km, data=data_div, family= Gamma (link = "identity"))
a3kq <- glm(alpha_FS ~ fc_3km + I(fc_3km^2), data=data_div, family= Gamma (link = "identity"))
## AIC-based model selection
(a.fs.aic <- AICctab(a0, a5k, a5kq,
a3k, a3kq,
base=T,weights=T, logLik=T))## logLik AICc dLogLik dAICc df weight
## a5k -19.3 47.6 2.8 0.0 3 0.523
## a3k -20.3 49.5 1.9 1.9 3 0.199
## a0 -22.1 49.6 0.0 2.0 2 0.195
## a5kq -19.0 51.8 3.1 4.2 4 0.064
## a3kq -20.3 54.2 1.9 6.6 4 0.019
Summary top model
summary(a5k)##
## Call:
## glm(formula = alpha_FS ~ fc_5km, family = Gamma(link = "identity"),
## data = data_div)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -0.43224 -0.11774 -0.03674 0.08987 0.35754
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 7.42155 0.94502 7.853 1.38e-05 ***
## fc_5km -0.07200 0.02815 -2.558 0.0285 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for Gamma family taken to be 0.06034959)
##
## Null deviance: 0.97549 on 11 degrees of freedom
## Residual deviance: 0.61269 on 10 degrees of freedom
## AIC: 44.608
##
## Number of Fisher Scoring iterations: 4
a.fs <- a5k# rename top model to save and plot belowInspect residuals:
sim <- simulateResiduals(a5k)
plot(sim)testDispersion(sim)##
## DHARMa nonparametric dispersion test via sd of residuals fitted vs.
## simulated
##
## data: simulationOutput
## dispersion = 1.1043, p-value = 0.656
## alternative hypothesis: two.sided
hnp::hnp(a5k)## Gamma model
boot::glm.diag.plots(a5k)Prediction:
Beta diversity additive
Modelling
### LM Beta for regular model selection
b0 <- glm(betad_FS~ 1, data=data_div, family= Gamma (link = "identity"))
b5k <- glm(betad_FS ~ fc_5km, data=data_div, family= Gamma (link = "identity"))
b5kq <- glm(betad_FS ~ fc_5km + I(fc_5km^2), data=data_div, family= Gamma (link = "identity"))
b3k <- glm(betad_FS ~ fc_3km, data=data_div, family= Gamma (link = "identity"))
b3kq <- glm(betad_FS ~ fc_3km + I(fc_3km^2), data=data_div, family= Gamma (link = "identity"))
## AIC-based model selection
(b.fs.aic <- AICctab(b0, b5k, b5kq,
b3k, b3kq,
base=T,weights=T, logLik=T))## logLik AICc dLogLik dAICc df weight
## b5k -22.1 53.3 5.1 0.0 3 0.690
## b5kq -21.4 56.6 5.8 3.3 4 0.133
## b3k -23.9 56.8 3.3 3.6 3 0.116
## b3kq -22.7 59.2 4.5 5.9 4 0.036
## b0 -27.2 59.8 0.0 6.6 2 0.026
Summary top model
summary(b5k)##
## Call:
## glm(formula = betad_FS ~ fc_5km, family = Gamma(link = "identity"),
## data = data_div)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -0.24293 -0.12664 -0.05680 0.08772 0.34333
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 13.29031 1.24634 10.663 8.8e-07 ***
## fc_5km -0.14349 0.03623 -3.961 0.00268 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for Gamma family taken to be 0.03508083)
##
## Null deviance: 0.76409 on 11 degrees of freedom
## Residual deviance: 0.32769 on 10 degrees of freedom
## AIC: 50.257
##
## Number of Fisher Scoring iterations: 4
b.fs <- b5k # rename top model to save and plot belowInspect residuals:
sim <- simulateResiduals(b5k)
plot(sim)testDispersion(sim) # ok##
## DHARMa nonparametric dispersion test via sd of residuals fitted vs.
## simulated
##
## data: simulationOutput
## dispersion = 1.1585, p-value = 0.608
## alternative hypothesis: two.sided
hnp::hnp(b5k) ## for glm Gamma or neg binom (maybe also others)## Gamma model
boot::glm.diag.plots(b5k) ## for glm Gamma or neg binomPrediction
bRaup-Crick (presence/absence-based)
Modelling
### LM gamma for regular model selection
rc0 <- lm(brc.fsNA ~ 1, data=data_div)
rc5k <- lm(brc.fsNA ~ fc_5km, data=data_div)
rc5kq <- lm(brc.fsNA ~ fc_5km + I(fc_5km^2), data=data_div)
rc3k <- lm(brc.fsNA ~ fc_3km, data=data_div)
rc3kq <- lm(brc.fsNA ~ fc_3km + I(fc_3km^2), data=data_div)
## AIC-based model selection
(rc.fs.aic <- AICctab(rc0, rc5k, rc5kq,
rc3k, rc3kq,
base=T,weights=T, logLik=T))## logLik AICc dLogLik dAICc df weight
## rc0 5.4 -5.5 0.0 0.0 2 0.482
## rc3kq 8.7 -3.7 3.3 1.8 4 0.197
## rc5k 5.9 -2.8 0.5 2.7 3 0.126
## rc5kq 8.1 -2.5 2.7 3.0 4 0.108
## rc3k 5.5 -2.1 0.1 3.4 3 0.087
Summary of 2nd best model, as the top is of absence of effect (~1)
summary(rc3kq)##
## Call:
## lm(formula = brc.fsNA ~ fc_3km + I(fc_3km^2), data = data_div)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.13558 -0.09064 -0.03247 0.07741 0.23095
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -0.6952262 0.2224103 -3.126 0.0122 *
## fc_3km 0.0409921 0.0173432 2.364 0.0424 *
## I(fc_3km^2) -0.0007333 0.0002926 -2.507 0.0335 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.1353 on 9 degrees of freedom
## Multiple R-squared: 0.4227, Adjusted R-squared: 0.2944
## F-statistic: 3.295 on 2 and 9 DF, p-value: 0.08439
rc1.fs <- rc0# rename top model to save and plot belowInspect residuals:
sim <- simulateResiduals(rc3kq)
plot(sim)testDispersion(sim) # ok##
## DHARMa nonparametric dispersion test via sd of residuals fitted vs.
## simulated
##
## data: simulationOutput
## dispersion = 0.82544, p-value = 0.816
## alternative hypothesis: two.sided
hnp::hnp(rc3kq)## Gaussian model (lm object)
Prediction:
Total abundance
Modelling
### LM "gamma" abundanec for regular model selection
gab0 <- glm.nb(ab.land_FS ~ 1, data=data_div)
gab5k <- glm.nb(ab.land_FS ~ fc_5km, data=data_div)
gab5kq <- glm.nb(ab.land_FS ~ fc_5km + I(fc_5km^2), data=data_div)
gab3k <- glm.nb(ab.land_FS ~ fc_3km, data=data_div)
gab3kq <- glm.nb(ab.land_FS ~ fc_3km + I(fc_3km^2), data=data_div)
## AIC-based model selection
(gab.fs.aic <- AICctab(gab0, gab5k, gab5kq,
gab3k, gab3kq,
base=T,weights=T, logLik=T))## logLik AICc dLogLik dAICc df weight
## gab5k -70.5 150.1 2.1 0.0 3 0.399
## gab0 -72.6 150.5 0.0 0.5 2 0.318
## gab3k -71.2 151.3 1.4 1.2 3 0.217
## gab5kq -70.4 154.5 2.2 4.5 4 0.043
## gab3kq -71.0 155.8 1.6 5.7 4 0.023
Summary 2nd best model, as top model was the one of absence of effect (~1):
summary(gab5k)##
## Call:
## glm.nb(formula = ab.land_FS ~ fc_5km, data = data_div, init.theta = 4.537178159,
## link = log)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -1.6496 -0.9873 -0.2873 0.6438 1.8333
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 5.93023 0.31931 18.572 <2e-16 ***
## fc_5km -0.02485 0.01110 -2.239 0.0251 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for Negative Binomial(4.5372) family taken to be 1)
##
## Null deviance: 17.268 on 11 degrees of freedom
## Residual deviance: 12.411 on 10 degrees of freedom
## AIC: 147.08
##
## Number of Fisher Scoring iterations: 1
##
##
## Theta: 4.54
## Std. Err.: 1.83
##
## 2 x log-likelihood: -141.082
gab.fs <- gab5k# rename top model to save and plot belowInspect residuals:
sim <- simulateResiduals(gab5k)
plot(sim)testDispersion(sim)##
## DHARMa nonparametric dispersion test via sd of residuals fitted vs.
## simulated
##
## data: simulationOutput
## dispersion = 1.1614, p-value = 0.56
## alternative hypothesis: two.sided
hnp::hnp(gab5k)## Negative binomial model (using MASS package)
boot::glm.diag.plots(gab5k)# R2 based on Nakagawa et al. 2017: variação explicada pelos efeitos fixos dividida pela total
gab5k.r2 <- r2(gab5k)
gab3k.r2 <- r2(gab3k)
## build AIC table do add model selection results:
aic.tab <- rbind(aic.tab, c("Total abundance FS",
"FC at 5km",
round(gab.fs.aic$AICc[1],2),
round(gab.fs.aic$dAICc[1],1),
gab.fs.aic$df[1],
round(gab.fs.aic$weight[1],2),
round(gab5k.r2$R2_Nagelkerke,2)),
c(" ", "~ 1 (NULL)",
round(gab.fs.aic$AICc[2],2),
round(gab.fs.aic$dAICc[2],1),
gab.fs.aic$df[2],
round(gab.fs.aic$weight[2],2),
"-"),
c(" ", "FC at 3km",
round(gab.fs.aic$AICc[3],2),
round(gab.fs.aic$dAICc[3],1),
gab.fs.aic$df[3],
round(gab.fs.aic$weight[3],2),
round(gab3k.r2$R2_Nagelkerke,2)))Prediction:
bRaup-Crick (abundance-based)
Modelling:
### LM gamma for regular model selection
rc0 <- lm(brc.abund.fsNA ~ 1, data=data_div)
rc5k <- lm(brc.abund.fsNA ~ fc_5km, data=data_div)
rc5kq <- lm(brc.abund.fsNA ~ fc_5km + I(fc_5km^2), data=data_div)
rc3k <- lm(brc.abund.fsNA ~ fc_3km, data=data_div)
rc3kq <- lm(brc.abund.fsNA ~ fc_3km + I(fc_3km^2), data=data_div)
## AIC-based model selection
(rc.fs.aic <- AICctab(rc0, rc5k, rc5kq,
rc3k, rc3kq,
base=T,weights=T, logLik=T))## logLik AICc dLogLik dAICc df weight
## rc5k 0.2 8.6 3.8 0.0 3 0.596
## rc3k -0.7 10.4 2.9 1.9 3 0.235
## rc0 -3.6 12.5 0.0 3.9 2 0.085
## rc5kq 0.3 13.1 3.9 4.5 4 0.062
## rc3kq -0.7 15.1 2.9 6.6 4 0.022
Summary top model:
summary(rc5k)##
## Call:
## lm(formula = brc.abund.fsNA ~ fc_5km, data = data_div)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.39341 -0.17488 -0.01206 0.15318 0.50889
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.90123 0.17499 5.150 0.000431 ***
## fc_5km -0.01796 0.00606 -2.964 0.014189 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.2602 on 10 degrees of freedom
## Multiple R-squared: 0.4677, Adjusted R-squared: 0.4144
## F-statistic: 8.786 on 1 and 10 DF, p-value: 0.01419
rc2.fs <- rc5k# rename top model to save and plot belowInspect residuals:
sim <- simulateResiduals(rc5k)
plot(sim)testDispersion(sim) # ok##
## DHARMa nonparametric dispersion test via sd of residuals fitted vs.
## simulated
##
## data: simulationOutput
## dispersion = 0.91715, p-value = 0.984
## alternative hypothesis: two.sided
hnp::hnp(rc5k)## Gaussian model (lm object)
Prediction
NON-FOREST SPECIALISTS
Gamma diversity
Modelling:
### LM alpha for regular model selection
g0 <- glm(gamma_NFS ~ 1, data=data_div, family= Gamma (link = "identity"))
g5k <- glm(gamma_NFS ~ fc_5km, data=data_div, family= Gamma (link = "identity"))
g5kq <- glm(gamma_NFS ~ fc_5km + I(fc_5km^2), data=data_div, family= Gamma (link = "identity"))
g3k <- glm(gamma_NFS ~ fc_3km, data=data_div, family= Gamma (link = "identity"))
g3kq <- glm(gamma_NFS ~ fc_3km + I(fc_3km^2), data=data_div, family= Gamma (link = "identity"))
## AIC-based model selection
(g.nfs.aic <- AICctab(g0,
g5k, g5kq,
g3k, g3kq,
base=T,weights=T, logLik=T))## logLik AICc dLogLik dAICc df weight
## g3k -27.3 63.5 7.6 0.0 3 0.7911
## g5k -29.2 67.5 5.7 4.0 3 0.1095
## g3kq -27.1 68.0 7.8 4.4 4 0.0864
## g5kq -29.2 72.2 5.7 8.6 4 0.0105
## g0 -34.9 75.2 0.0 11.6 2 0.0024
Summary top model:
summary(g3k)##
## Call:
## glm(formula = gamma_NFS ~ fc_3km, family = Gamma(link = "identity"),
## data = data_div)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -0.38392 -0.22774 -0.05573 0.15121 0.43058
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 18.37185 2.13069 8.622 6.07e-06 ***
## fc_3km -0.30280 0.05227 -5.793 0.000175 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for Gamma family taken to be 0.08050867)
##
## Null deviance: 2.71637 on 11 degrees of freedom
## Residual deviance: 0.78116 on 10 degrees of freedom
## AIC: 60.54
##
## Number of Fisher Scoring iterations: 4
g.nfs <- g3k # rename top model to save and plot belowInspect residuals:
sim <- simulateResiduals(g3k)
plot(sim)testDispersion(sim) # a lil under##
## DHARMa nonparametric dispersion test via sd of residuals fitted vs.
## simulated
##
## data: simulationOutput
## dispersion = 0.90839, p-value = 0.96
## alternative hypothesis: two.sided
hnp::hnp(g3k) ## for glm poisson or neg binom (maybe also others)## Gamma model
boot::glm.diag.plots(g3k) ## for glm poisson or neg binomPrediction:
Alpha diversity
Modelling:
### LM alpha for regular model selection
a0 <- glm(alpha_NFS ~ 1, data=data_div, family= Gamma (link = "identity"))
a5k <- glm(alpha_NFS ~ fc_5km, data=data_div, family= Gamma (link = "identity"))
a5kq <- glm(alpha_NFS ~ fc_5km + I(fc_5km^2), data=data_div, family= Gamma (link = "identity"))
a3k <- glm(alpha_NFS ~ fc_3km, data=data_div, family= Gamma (link = "identity"))
a3kq <- glm(alpha_NFS ~ fc_3km + I(fc_3km^2), data=data_div, family= Gamma (link = "identity"))
## AIC-based model selection
(a.nfs.aic <- AICctab(a0, a5k, a5kq,
a3k, a3kq,
base=T,weights=T, logLik=T))## logLik AICc dLogLik dAICc df weight
## a3k -16.7 42.3 5.7 0.0 3 0.490
## a5k -16.9 42.9 5.4 0.6 3 0.371
## a5kq -16.2 46.2 6.1 3.9 4 0.071
## a3kq -16.5 46.6 5.9 4.3 4 0.057
## a0 -22.3 50.0 0.0 7.7 2 0.011
Summary top model:
summary(a3k)##
## Call:
## glm(formula = alpha_NFS ~ fc_3km, family = Gamma(link = "identity"),
## data = data_div)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -0.68238 -0.15994 0.00108 0.11140 0.55352
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 5.82215 0.83950 6.935 4.02e-05 ***
## fc_3km -0.09207 0.02113 -4.357 0.00143 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for Gamma family taken to be 0.1163241)
##
## Null deviance: 3.1391 on 11 degrees of freedom
## Residual deviance: 1.2537 on 10 degrees of freedom
## AIC: 39.336
##
## Number of Fisher Scoring iterations: 6
a.nfs <- a3k # rename top model to save and plot belowInspect residuals:
sim <- simulateResiduals(a3k)
plot(sim)testDispersion(sim) # a lil under##
## DHARMa nonparametric dispersion test via sd of residuals fitted vs.
## simulated
##
## data: simulationOutput
## dispersion = 0.93258, p-value = 0.928
## alternative hypothesis: two.sided
hnp::hnp(a3k) ## for glm poisson or neg binom (maybe also others)## Gamma model
boot::glm.diag.plots(a3k) ## for glm poisson or neg binomPrediction:
Beta diversity additive
Modelling
### LM Beta for regular model selection
b0 <- glm(betad_NFS~ 1, data=data_div, family= Gamma (link = "identity"))
b5k <- glm(betad_NFS ~ fc_5km, data=data_div, family= Gamma (link = "identity"))
b5kq <- glm(betad_NFS ~ fc_5km + I(fc_5km^2), data=data_div, family= Gamma (link = "identity"))
b3k <- glm(betad_NFS ~ fc_3km, data=data_div, family= Gamma (link = "identity"))
b3kq <- glm(betad_NFS ~ fc_3km + I(fc_3km^2), data=data_div, family= Gamma (link = "identity"))
## AIC-based model selection
(b.nfs.aic <- AICctab(b0, b5k, b5kq,
b3k, b3kq,
base=T,weights=T, logLik=T))## logLik AICc dLogLik dAICc df weight
## b3k -25.0 58.9 5.9 0.0 3 0.753
## b5k -26.8 62.7 4.0 3.7 3 0.116
## b3kq -24.6 62.9 6.3 3.9 4 0.105
## b5kq -26.6 67.0 4.2 8.1 4 0.013
## b0 -30.9 67.1 0.0 8.2 2 0.013
Summary top model:
summary(b3k)##
## Call:
## glm(formula = betad_NFS ~ fc_3km, family = Gamma(link = "identity"),
## data = data_div)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -0.71147 -0.22545 -0.03265 0.18829 0.45483
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 12.59989 1.68741 7.467 2.14e-05 ***
## fc_3km -0.21226 0.04078 -5.205 0.000399 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for Gamma family taken to be 0.1116167)
##
## Null deviance: 3.2119 on 11 degrees of freedom
## Residual deviance: 1.2312 on 10 degrees of freedom
## AIC: 55.933
##
## Number of Fisher Scoring iterations: 4
b.nfs <- b3k # rename top model to save and plot belowInspect residuals:
sim <- simulateResiduals(b3k)
plot(sim)testDispersion(sim)##
## DHARMa nonparametric dispersion test via sd of residuals fitted vs.
## simulated
##
## data: simulationOutput
## dispersion = 0.65877, p-value = 0.68
## alternative hypothesis: two.sided
hnp::hnp(b3k) ## for glm poisson or neg binom (maybe also others)## Gamma model
boot::glm.diag.plots(b3k) ## for glm poisson or neg binomPrediction:
bRaup-Crick (presence/absence-based)
Modelling
### LM gamma for regular model selection
rc0 <- lm(brc.nfsNA ~ 1, data=data_div)
rc5k <- lm(brc.nfsNA ~ fc_5km, data=data_div)
rc5kq <- lm(brc.nfsNA ~ fc_5km + I(fc_5km^2), data=data_div)
rc3k <- lm(brc.nfsNA ~ fc_3km, data=data_div)
rc3kq <- lm(brc.nfsNA ~ fc_3km + I(fc_3km^2), data=data_div)
## AIC-based model selection
(rc.nfs.aic <- AICctab(rc0, rc5k, rc5kq,
rc3k, rc3kq,
base=T,weights=T, logLik=T))## logLik AICc dLogLik dAICc df weight
## rc3k -0.5 10.0 2.3 0.0 3 0.393
## rc0 -2.8 11.0 0.0 1.0 2 0.244
## rc5k -1.1 11.2 1.7 1.2 3 0.220
## rc3kq 0.6 12.5 3.4 2.4 4 0.116
## rc5kq -0.8 15.4 2.0 5.3 4 0.027
Summary of 2nd best model, 1st wqs of absence of effect (~1):
summary(rc3k)##
## Call:
## lm(formula = brc.nfsNA ~ fc_3km, data = data_div)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.32362 -0.21530 -0.05601 0.23870 0.39743
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.32268 0.18990 1.699 0.1201
## fc_3km -0.01402 0.00647 -2.167 0.0555 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.2766 on 10 degrees of freedom
## Multiple R-squared: 0.3195, Adjusted R-squared: 0.2515
## F-statistic: 4.696 on 1 and 10 DF, p-value: 0.05545
rc1.nfs <- rc3k # rename top model to save and plot belowInspect residuals:
sim <- simulateResiduals(rc3k)
plot(sim)testDispersion(sim) # a lil under##
## DHARMa nonparametric dispersion test via sd of residuals fitted vs.
## simulated
##
## data: simulationOutput
## dispersion = 0.91715, p-value = 0.984
## alternative hypothesis: two.sided
hnp::hnp(rc3k)## Gaussian model (lm object)
Prediction:
Total abundance
Modelling:
### LM "gamma" abundanec for regular model selection
gab0 <- glm.nb(ab.land_NFS ~ 1, data=data_div)
gab5k <- glm.nb(ab.land_NFS ~ fc_5km, data=data_div)
gab5kq <- glm.nb(ab.land_NFS ~ fc_5km + I(fc_5km^2), data=data_div)
gab3k <- glm.nb(ab.land_NFS ~ fc_3km, data=data_div)
gab3kq <- glm.nb(ab.land_NFS ~ fc_3km + I(fc_3km^2), data=data_div)
## AIC-based model selection
(gab.nfs.aic <- AICctab(gab0, gab5k, gab5kq,
gab3k, gab3kq,
base=T,weights=T, logLik=T))## logLik AICc dLogLik dAICc df weight
## gab5k -63.9 136.9 5.3 0.0 3 0.404
## gab3k -64.3 137.6 4.9 0.7 3 0.284
## gab5kq -62.1 137.8 7.2 0.9 4 0.252
## gab3kq -63.7 141.2 5.5 4.3 4 0.047
## gab0 -69.2 143.8 0.0 6.9 2 0.013
Summary top model:
summary(gab5k)##
## Call:
## glm.nb(formula = ab.land_NFS ~ fc_5km, data = data_div, init.theta = 4.164123767,
## link = log)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -2.8516 -0.4183 0.2175 0.4846 0.8634
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 5.93704 0.33661 17.638 < 2e-16 ***
## fc_5km -0.04838 0.01182 -4.094 4.24e-05 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for Negative Binomial(4.1641) family taken to be 1)
##
## Null deviance: 29.266 on 11 degrees of freedom
## Residual deviance: 12.765 on 10 degrees of freedom
## AIC: 133.88
##
## Number of Fisher Scoring iterations: 1
##
##
## Theta: 4.16
## Std. Err.: 1.77
##
## 2 x log-likelihood: -127.88
gab.nfs <- gab5k # rename top model to save and plot belowInspect residuals:
sim <- simulateResiduals(gab5k)
plot(sim)testDispersion(sim) # a lil under##
## DHARMa nonparametric dispersion test via sd of residuals fitted vs.
## simulated
##
## data: simulationOutput
## dispersion = 0.45825, p-value = 0.288
## alternative hypothesis: two.sided
hnp::hnp(gab5k)## Negative binomial model (using MASS package)
boot::glm.diag.plots(gab5k)# R2 based on Nakagawa et al. 2017: variação explicada pelos efeitos fixos dividida pela total
gab5k.r2 <- r2(gab5k)
gab5kq.r2 <- r2(gab5kq)
gab3k.r2 <- r2(gab3k)
## build AIC table do add model selection results:
aic.tab <- rbind(aic.tab, c("Total abundance NFS",
"FC at 5km",
round(gab.nfs.aic$AICc[1],2),
round(gab.nfs.aic$dAICc[1],1),
gab.nfs.aic$df[1],
round(gab.nfs.aic$weight[1],2),
round(gab5k.r2$R2_Nagelkerke,2)),
c(" ", "FC at 3km",
round(gab.nfs.aic$AICc[2],2),
round(gab.nfs.aic$dAICc[2],1),
gab.nfs.aic$df[2],
round(gab.nfs.aic$weight[2],2),
round(gab3k.r2$R2_Nagelkerke,2)),
c(" ", "Non-linear FC at 5km",
round(gab.nfs.aic$AICc[3],2),
round(gab.nfs.aic$dAICc[3],1),
gab.nfs.aic$df[3],
round(gab.nfs.aic$weight[3],2),
round(gab5kq.r2$R2_Nagelkerke,2)
))Prediction:
bRaup-Crick (abundance-based)
Modelling
### LM gamma for regular model selection
rc0 <- lm(brc.abund.nfsNA ~ 1, data=data_div)
rc5k <- lm(brc.abund.nfsNA ~ fc_5km, data=data_div)
rc5kq <- lm(brc.abund.nfsNA ~ fc_5km + I(fc_5km^2), data=data_div)
rc3k <- lm(brc.abund.nfsNA ~ fc_3km, data=data_div)
rc3kq <- lm(brc.abund.nfsNA ~ fc_3km + I(fc_3km^2), data=data_div)
## AIC-based model selection
(rc.nfs.aic <- AICctab(rc0, rc5k, rc5kq,
rc3k, rc3kq,
base=T,weights=T, logLik=T))## logLik AICc dLogLik dAICc df weight
## rc3k -7.5 24.0 1.9 0.0 3 0.377
## rc0 -9.4 24.2 0.0 0.2 2 0.341
## rc5k -8.0 25.1 1.4 1.1 3 0.219
## rc3kq -7.4 28.5 2.1 4.5 4 0.041
## rc5kq -8.0 29.7 1.5 5.7 4 0.022
Summary 2nd best model, 1st is of absence of effect (~1):
summary(rc3k)##
## Call:
## lm(formula = brc.abund.nfsNA ~ fc_3km, data = data_div)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1.25085 -0.18926 0.04247 0.27112 0.52852
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.84365 0.33999 2.481 0.0325 *
## fc_3km -0.02260 0.01158 -1.951 0.0797 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.4951 on 10 degrees of freedom
## Multiple R-squared: 0.2756, Adjusted R-squared: 0.2032
## F-statistic: 3.805 on 1 and 10 DF, p-value: 0.07965
rc2.nfs <- rc3k # rename top model to save and plot belowInspect residuals:
sim <- simulateResiduals(rc3k)
plot(sim)testDispersion(sim) # a lil under##
## DHARMa nonparametric dispersion test via sd of residuals fitted vs.
## simulated
##
## data: simulationOutput
## dispersion = 0.91715, p-value = 0.984
## alternative hypothesis: two.sided
hnp::hnp(rc3k)## Gaussian model (lm object)
Prediction:
AICc-model selection table
| Diversity | Model | AICc | dAICc | df | weight | r2 |
|---|---|---|---|---|---|---|
| Gamma FS | FC at 5km | 64.84 | 0 | 3 | 0.73 | 0.55 |
| Alpha FS | FC at 5km | 47.61 | 0 | 3 | 0.52 | 0.38 |
| FC at 3km | 49.54 | 1.9 | 3 | 0.2 | 0.27 | |
| ~ 1 (NULL) | 49.58 | 2 | 2 | 0.19 |
|
|
| Beta FS | FC at 5km | 53.26 | 0 | 3 | 0.69 | 0.58 |
| RCbeta (P/A-based) FS | ~ 1 (NULL) | -5.48 | 0 | 2 | 0.48 |
|
| Non-linear FC at 3km | -3.7 | 1.8 | 4 | 0.2 | 0.29 | |
| Total abundance FS | FC at 5km | 150.08 | 0 | 3 | 0.4 | 0.44 |
| ~ 1 (NULL) | 150.54 | 0.5 | 2 | 0.32 |
|
|
| FC at 3km | 151.31 | 1.2 | 3 | 0.22 | 0.32 | |
| RCbeta (abund-based) FS | FC at 5km | 8.56 | 0 | 3 | 0.6 | 0.41 |
| FC at 3km | 10.41 | 1.9 | 3 | 0.24 | 0.32 | |
| Gamma NFS | FC at 3km | 63.54 | 0 | 3 | 0.79 | 0.74 |
| Alpha NFS | FC at 3km | 42.34 | 0 | 3 | 0.49 | 0.63 |
| FC at 5km | 42.89 | 0.6 | 3 | 0.37 | 0.61 | |
| Beta NFS | FC at 3km | 58.93 | 0 | 3 | 0.75 | 0.65 |
| RCbeta (P/A-based) NFS | FC at 3km | 10.02 | 0 | 3 | 0.39 | 0.25 |
| ~ 1 (NULL) | 10.97 | 1 | 2 | 0.24 |
|
|
| FC at 5km | 11.18 | 1.2 | 3 | 0.22 | 0.18 | |
| Total abundance NFS | FC at 5km | 136.88 | 0 | 3 | 0.4 | 0.82 |
| FC at 3km | 137.58 | 0.7 | 3 | 0.28 | 0.79 | |
| Non-linear FC at 5km | 137.82 | 0.9 | 4 | 0.25 | 0.93 | |
| RCbeta (abund-based) NFS | FC at 3km | 24 | 0 | 3 | 0.38 | 0.2 |
| ~ 1 (NULL) | 24.2 | 0.2 | 2 | 0.34 |
|
|
| FC at 5km | 25.09 | 1.09 | 3 | 0.22 | 0.13 |
FIGURES
Gamma hab
### Gamma ####
# r2(g.fs) # 0.55
## prediction g2
new <- data.frame(fc_5km = seq(min(data_div$fc_5km),max(data_div$fc_5km), length.out=50))
preds <- predict(g.fs, type="response", newdata=new, se.fit=T)
new$pred <- preds$fit
new$ic.up <- preds$fit + preds$se.fit*1.96
new$ic.low <- preds$fit - preds$se.fit*1.96
# r2(g.nfs) # 0.73
## prediction g1
new2 <- data.frame(fc_3km = seq(min(data_div$fc_3km),max(data_div$fc_3km), length.out=50))
preds2 <- predict(g.nfs, type="response", newdata=new2, se.fit=T)
new2$pred <- preds2$fit
new2$ic.up <- preds2$fit + preds2$se.fit*1.96
new2$ic.low <- preds2$fit - preds2$se.fit*1.96
## Plot gamma richness ~ fc %
(g.hab <- ggplot(data= new, aes(x= fc_5km, y = pred)) +
geom_line(data= new, aes(x=fc_5km, y=pred), col= "#0072B2") +
geom_ribbon(data= new, aes(x=fc_5km , ymax=ic.up, ymin=ic.low), alpha=0.1, fill= "#0072B2") +
geom_point(data= data_div, aes(x=fc_5km, y=gamma_FS), col= "#0072B2") +
geom_ribbon(data= new2, aes(x=fc_3km , ymax=ic.up, ymin=ic.low), alpha=0.1, fill= "#D55E00") +
geom_line(data= new2, aes(x=fc_3km, y=pred), col= "#D55E00") +
geom_point(data= data_div, aes(x=fc_3km, y=gamma_NFS), col= "#D55E00") +
xlab(" ") + ylab("Gamma") +
theme(axis.text=element_text(size=15),axis.title=element_text(size=14,face="bold"),legend.position = 'none') +
ylim(0,max(data_div[2:7])+0.5) +
# scale_x_reverse() +
annotate("text", y=20.5, x=38,colour= "#0072B2",
label=expression(paste(~ R^2 ~ "= 0.55"))) +
annotate("text", y=19, x=40, colour= "#D55E00",
label=expression(paste(~ R^2 ~ "= 0.73"))) + labs(tag="a"))Alpha hab
# r2(a.fs) # 0.38
## prediction a2
new <- data.frame(fc_5km = seq(min(data_div$fc_5km),max(data_div$fc_5km), length.out=50))
preds <- predict(a.fs, type="response", newdata=new, se.fit=T)
new$pred <- preds$fit
new$ic.up <- preds$fit + preds$se.fit*1.96
new$ic.low <- preds$fit - preds$se.fit*1.96
# r2(a.nfs) # 0.63
## prediction a1
new2 <- data.frame(fc_3km = seq(min(data_div$fc_3km),max(data_div$fc_3km), length.out=50))
preds <- predict(a.nfs, type="response", newdata=new2, se.fit=T)
new2$pred <- preds$fit
new2$ic.up <- preds$fit + preds$se.fit*1.96
new2$ic.low <- preds$fit - preds$se.fit*1.96
## Plot Alpha richness ~ fc_3km
(a.hab <- ggplot(data= new, aes(x=fc_5km, y=pred)) +
geom_line(data= new, aes(x=fc_5km, y=pred), col= "#0072B2",
alpha= .5, linetype= "dashed") +
# geom_ribbon(data= new, aes(x=fc_5km , ymax=ic.up, ymin=ic.low), alpha=0.1, fill= "#0072B2") +
geom_point(data=data_div, aes(x=fc_5km, y=alpha_FS), col= "#0072B2") +
geom_line(data= new2, aes(x=fc_3km, y=pred), col= "#D55E00") +
geom_ribbon(data= new2, aes(x=fc_3km , ymax=ic.up, ymin=ic.low), alpha=0.1, fill= "#D55E00") +
geom_point(data=data_div, aes(x=fc_3km, y=alpha_NFS), col= "#D55E00") +
xlab("") + ylab("Alpha") +
scale_color_manual(breaks=c("FS", "NFS"),
values=c('FS'='#0072B2', 'NFS'='#D55E00')) +
theme(axis.text=element_text(size=15),
axis.title=element_text(size=14,face="bold"),
legend.position = 'bottom') +
ylim(0,max(data_div[2:7])+0.5) +
# scale_x_reverse() +
annotate("text", y=13, x=38,colour= "#0072B2",
label=expression(paste(~ R^2 ~ "= 0.38")))+
annotate("text", y=11, x=40, colour= "#D55E00",
label=expression(paste(~ R^2 ~ "= 0.63"))) +
labs(tag="b"))Beta hab
### Beta add ####
# r2(b.fs) # 0.58
## prediction b1
new <- data.frame(fc_5km = seq(min(data_div$fc_5km),max(data_div$fc_5km), length.out=50))
preds <- predict(b.fs, type="response", newdata=new, se.fit=T)
new$pred <- preds$fit
new$ic.up <- preds$fit + preds$se.fit*1.96
new$ic.low <- preds$fit - preds$se.fit*1.96
# r2(b.nfs) # 0.65
## prediction b1
new2 <- data.frame(fc_3km = seq(min(data_div$fc_3km),max(data_div$fc_3km), length.out=50))
preds <- predict(b.nfs, type="response", newdata=new2, se.fit=T)
new2$pred <- preds$fit
new2$ic.up <- preds$fit + preds$se.fit*1.96
new2$ic.low <- preds$fit - preds$se.fit*1.96
## Plot Beta Diversity ~ fc %
(b.hab <- ggplot(data= new, aes(x=fc_5km, y=pred)) +
geom_line(data= new, aes(x=fc_5km, y=pred), col= "#0072B2") +
geom_ribbon(data= new, aes(x=fc_5km , ymax=ic.up, ymin=ic.low), alpha=0.1, fill= "#0072B2") +
geom_point(data=data_div, aes(x=fc_5km, y=betad_FS), col= "#0072B2") +
geom_line(data= new2, aes(x= fc_3km, y= pred),
col= "#D55E00") +
geom_ribbon(data= new2, aes(x=fc_3km , ymax=ic.up, ymin=ic.low), alpha=0.1, fill= "#D55E00") +
geom_point(data=data_div, aes(x=fc_3km, y=betad_NFS), col= "#D55E00") +
xlab(" ") + ylab("Beta") +
theme(axis.text=element_text(size=15),axis.title=element_text(size=14,face="bold"),legend.position = 'none') +
ylim(0,max(data_div[2:7])+0.5) +
# scale_x_reverse() +
annotate("text", y=18, x=40,colour= "#0072B2",
label=expression(paste(~ R^2 ~ "= 0.58")))+
annotate("text", y=16, x=38, colour= "#D55E00",
label=expression(paste(~ R^2 ~ "= 0.65"))) +
labs(tag="c"))g.hab + a.hab + b.habggsave("figures/plot.hab.jpeg",units="cm", width=20, height=8, dpi=300)Total abundance FS & NFS
## FS
# r2(gab.fs) # 0.44
## prediction gab1
new <- data.frame(fc_5km = seq(min(data_div$fc_5km),max(data_div$fc_5km), length.out=50))
preds <- predict(gab.fs, type="response", newdata=new, se.fit=T)
new$pred <- preds$fit
new$ic.up <- preds$fit + preds$se.fit*1.96
new$ic.low <- preds$fit - preds$se.fit*1.96
## NFS
# r2(gab.nfs) # 0.82
## prediction gab2
new2 <- new
preds <- predict(gab.nfs, type="response", newdata=new2, se.fit=T)
new2$pred <- preds$fit
new2$ic.up <- preds$fit + preds$se.fit*1.96
new2$ic.low <- preds$fit - preds$se.fit*1.96
## Plot gamma richness ~ forest cover 5km
(gab.hab <- ggplot(new, aes(x=fc_5km, y=pred)) +
geom_line(linetype= "dashed", col= "#0072B2", aes(alpha= 0.1)) +
# geom_ribbon(aes(x=fc_5km , ymax=ic.up, ymin=ic.low), alpha=0.1, fill= "#0072B2") +
geom_point(data=data_div, aes(x=fc_5km, y=ab.land_FS), col= "#0072B2") +
geom_line(data= new2, aes(x=fc_5km, y=pred), col= "#D55E00") +
geom_ribbon(data= new2, aes(x=fc_5km , ymax=ic.up, ymin=ic.low), alpha=0.2, fill= "#D55E00") +
geom_point(data= data_div, aes(x=fc_5km, y=ab.land_NFS), col= "#D55E00") +
scale_x_continuous(limits= c(10,50)) +
xlab(" ") + ylab("Total abundance") +
theme(axis.text=element_text(size=15),axis.title=element_text(size=14,face="bold"),legend.position = 'none') +
# scale_x_reverse() +
annotate("text", y=300, x=38, alpha= .5, colour= "#0072B2",
label=expression(paste(~ R^2 ~ "= 0.44"))) +
annotate("text", y=270, x=40, colour= "#D55E00",
label=expression(paste(~ R^2 ~ "= 0.82"))) +
labs(tag="d"))
### Beta RC FS & NFS
## FS
# ## NULL ## r2(rc1.fs)
## prediction rc2.2 and rc2
new <- data.frame(fc_3km = seq(min(data_div$fc_3km),max(data_div$fc_3km), length.out=30))
preds <- predict(rc1.fs, type="response", newdata=new, se.fit=T)
new$pred <- preds$fit
new$ic.up <- preds$fit + preds$se.fit*1.96
new$ic.low <- preds$fit - preds$se.fit*1.96
## NFS
# r2(rc1.nfs) # 0.25
## prediction rc2
new2 <- data.frame(fc_3km = seq(min(data_div$fc_3km),max(data_div$fc_3km), length.out=50))
preds <- predict(rc1.nfs, type="response", newdata=new2, se.fit=T)
new2$pred <- preds$fit
new2$ic.up <- preds$fit + preds$se.fit*1.96
new2$ic.low <- preds$fit - preds$se.fit*1.96
## Plot RCbeta abundance ~ forest cover at 3
(rc.hab <- ggplot(data= new, aes(x=fc_3km, y=pred)) +
geom_line(data= new,aes(x=fc_3km, y=pred),
col= "#0072B2", linetype= "dashed", alpha= .2) +
# geom_ribbon(data= new, aes(x=fc_3km , ymax=ic.up, ymin=ic.low), alpha=0.1, fill= "#0072B2") +
geom_point(data=data_div, aes(x=fc_3km, y=brc.fsNA), col= "#0072B2") +
geom_line(data= new2, aes(x=fc_3km, y=pred), col= "#D33E00", linetype= "dashed", alpha= .2) +
# geom_ribbon(data= new2, aes(x=fc_3km , ymax=ic.up, ymin=ic.low), alpha=0.1, fill= "#D33E00") +
geom_point(data=data_div, aes(x=fc_3km, y=brc.nfsNA), col= "#D33E00") +
xlab("Percent forest cover") + ylab(c(expression(bold("β"["RC"])))) +
theme(axis.text=element_text(size=14), axis.title=element_text(size=14,face="bold"),
legend.position = 'none') + ylim(-1,1) +
# scale_x_reverse() +
annotate("text", y=0.9, x=43, alpha= .5, colour= "#D33E00",
label=expression(paste(~ R^2 ~ "= 0.25"))) +
# annotate("text", y=0.7, x=43, alpha= .5, colour= "#0072B2",
# label=expression(paste(~ R^2 ~ "= 0.30"))) +
# annotate("text", y=0.9, x=43, colour= "#D33E00",
# label=expression(paste(beta, "(SE) = -0.06(0.03); " ~ R^2 ~ "= 0.37"))) +
# annotate("text", y=0.8, x=43, colour= "#0072B2",
# label=expression(paste(beta, "(SE) = -0.06(0.03); " ~ R^2 ~ "= 0.33"))) +
labs(tag="e"))
Beta RC abund FS & NFS
# r2(rc2.fs) # 0.41
## prediction rc2
new <- data.frame(fc_5km = seq(min(data_div$fc_5km),max(data_div$fc_5km), length.out=50))
preds <- predict(rc2.fs, type="response", newdata=new, se.fit=T)
new$pred <- preds$fit
new$ic.up <- preds$fit + preds$se.fit*1.96
new$ic.low <- preds$fit - preds$se.fit*1.96
## NFS
# r2(rc2.nfs) # 0.21
## prediction rc1
new2 <- data.frame(fc_3km = seq(min(data_div$fc_3km),max(data_div$fc_3km), length.out=50))
preds <- predict(rc2.nfs, type="response", newdata=new2, se.fit=T)
new2$pred <- preds$fit
new2$ic.up <- preds$fit + preds$se.fit*1.96
new2$ic.low <- preds$fit - preds$se.fit*1.96
## RCbeta abund ~ FC%
(rc.ab.hab <- ggplot(data= new, aes(x=fc_5km, y=pred)) +
geom_line(data= new,aes(x=fc_5km, y=pred), col= "#0072B2") +
geom_ribbon(data= new, aes(x=fc_5km , ymax=ic.up, ymin=ic.low), alpha=0.1, fill= "#0072B2") +
geom_point(data=data_div, aes(x=fc_5km, y=brc.abund.fsNA), col= "#0072B2") +
geom_line(data= new2, aes(x=fc_3km, y=pred), linetype= "dashed", col= "#D55E00", alpha=.2) +
# geom_ribbon(data= new2, aes(x=fc_3km , ymax=ic.up, ymin=ic.low), alpha=0.1, fill= "#D55E00") +
geom_point(data=data_div, aes(x=fc_3km, y=brc.abund.nfsNA), col= "#D55E00") +
xlab(" ") + ylab(c(expression(bold("β"["RC-abundance"])))) +
theme(axis.text=element_text(size=14), axis.title=element_text(size=14,face="bold"),legend.position = 'none') +
ylim(-1,1.1) +
# scale_x_reverse() +
annotate("text", y=0.9, x=45, colour= "#0072B2",
label=expression(paste(~ R^2 ~ "= 0.41"))) +
annotate("text", y=0.75, x=43, alpha= .5, colour= "#D55E00",
label=expression(paste( ~ R^2 ~ "= 0.21"))) +
labs(tag="f"))Final plot
(plot.hab <- (g.hab + a.hab + b.hab) /
(gab.hab + rc.hab + rc.ab.hab))ggsave("figures/figure2.jpeg",units="cm", width=28, height=15, dpi=300)